Tanzania
Rasmus Benestad
May 10, 2016
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
## Script to downscale temperature from Tanzania.
## Data from Habiba.
library(esd)Functions
Prepare the predictands
predictands <- function(param='tmin',incl.ghcn=FALSE) {
## Prepare the predictand data
#print("predictands")
fname <- paste(param,'.tz.rda',sep='')
path <- paste('Habiba/',param,sep='')
if (!file.exists(fname)) {
files <- list.files(path,pattern='rda',full.names=TRUE)
locs <- list.files(path,pattern='rda')
dse <- grep('dse.tm',files)
files <- files[-dse]; locs <- locs[-dse]
print(files)
locs <- substr(locs,1,nchar(locs)-4)
Y <- NULL
for (i in 1:length(files)) {
load(files[i])
y1 <- eval(parse(text=locs[i]))
if (is.null(Y)) Y <- y1 else
Y <- combine.stations(Y,y1)
}
save(file=fname,Y)
plot(Y)
} else load(fname)
#print('GHCN')
if (incl.ghcn) {
xstations <- switch(param,'tmin'='Habiba/tn.eafrica.rda',
'tmax'='Habiba/tx.eafrica.rda')
load(xstations)
if (param=='tmin') y <- tn.eafrica else y <- tx.eafrica
y <- subset(y,it=c(1950,2015))
Y <- combine(Y,y)
}
map(Y,FUN='mean',cex=-2)
invisible(Y)
}Function that applies the downscaling to a combination of predictands in terms of principal component analysis (PCA) results and predictors that involve large multi-model ensembles.
downscale <- function(Y,predictor,it='djf',param='t2m',FUN='mean',FUNX='mean',
period=c(1950,2015),plot=FALSE,rcp='rcp45',verbose=FALSE,
lon=c(0,100),lat=c(65,90),eofs=1:6,n=10,rel.cord=FALSE,select=NULL) {
print('downscale')
## Use a time and space window:
Y <- subset(Y,it=period)
Y <- subset(Y,is=list(lat=lat,lon=lon))
## Estimate seasonal means & weed out stations with little data
Y4 <- subset(as.4seasons(Y,FUN=FUN,nmin=60),it=it)
ok <- apply(coredata(Y4),1,nv)
Y4 <- subset(Y4,it=ok>0)
nok <- apply(coredata(Y4),2,nv)
Y4 <- subset(Y4,is=nok>15)
print(paste(round(100*sum(!is.finite(Y4))/length(!is.finite(Y4))),'% missing',sep=''))
if (plot) map(Y,FUN=FUN,cex=-2)
nmiss <- round(100*sum(!is.finite(Y4))/length(!is.finite(Y4)))
print(paste(nmiss,'% missing',sep=''))
## Fill missing data using PCA-based regression
Z <- pcafill(Y4)
pca <- PCA(Z,n=n)
if (plot) plot(pca)
## Downscale results
print('DSensemble')
dse.pca <- DSensemble(pca,predictor=predictor,FUNX=FUNX,verbose=verbose,
biascorrect=TRUE,rcp=rcp,eofs=eofs,select=select,
lon=lon,lat=lat,rel.cord=rel.cord,it=it)
attr(dse.pca,'N.missing') <- nmiss
invisible(dse.pca)
}Function that extracts results for a given year
distill <- function(x,it=2050) {
y <- subset(x,it=it)
attributes(y) <- NULL
y <- mean(y)
return(y)
}Function that grids the results
gridmap <- function(dse.results,param=NULL,
fast=TRUE,it=c(2070,2099),it0=c(1979,2012),verbose=FALSE,
breaks=seq(0,7,length.out=15),pal='warm',rev=FALSE) {
if (verbose) print('gridmap')
library(LatticeKrig)
if (is.null(param)) param <- varid(dse.results)
if (verbose) print ('as.station')
if (inherits(dse.results,'pca')) dse.station <- as.station(dse.results) else
dse.station <- dse.results
if (verbose) print('distill')
Z0 <- lapply(dse.station,distill,it=it0)
Z <- lapply(dse.station,distill,it=it)
lon <- unlist(lapply(dse.station,lon))
lat <- unlist(lapply(dse.station,lat))
alt <- unlist(lapply(dse.station,function(x) alt(attr(x,'station'))))
X <- attr(dse.station[[1]],'station')
z0 <- matrix(unlist(Z0),1,length(dse.station))
z <- matrix(unlist(Z),1,length(dse.station))
data(etopo5)
etopo5 <- subset(etopo5,is=list(lon=range(lon),lat=range(lat)))
etopo5[etopo5<=0] <- NA
if (verbose) print('grid')
## Set the grid to be the same as that of etopo5:
grid <- structure(list(x=lon(etopo5),y=lat(etopo5)),class='gridList')
#3 Flag dubplicated stations:
ok <- !(duplicated(lon) & duplicated(lat))
## Spread in the 90-percente interval changing
obj <- LatticeKrig( x=cbind(lon[ok],lat[ok]),
y=z[1,ok] - z0[1,ok],Z=alt[ok] )
## obj <- LatticeKrig( x=cbind(lon[ok],lat[ok]), y=z[2,ok],Z=alt[ok])
w <- predictSurface(obj, grid.list = grid,Z=etopo5)
w$z[is.na(etopo5)] <- NA
## Get rid of packages that have functions of same name:
detach("package:LatticeKrig")
detach("package:fields")
detach("package:spam")
detach("package:grid")
detach("package:maps")
## Convert the results from LatticeKrig to esd:
W <- w$z
attr(W,'variable') <- param
attr(W,'unit') <- 'degC'
attr(W,'longitude') <- w$x
attr(W,'latitude') <- w$y
class(W) <- class(etopo5)
## Make a projection that zooms in on the Barents region
#dev.new()
colbar <- list(breaks=round(breaks,2),rev=rev,pal=pal)
if (verbose) print(colbar)
if (is.null(rev)) rev <- switch(param,'t2m'=FALSE,'pr'=TRUE)
Wx <- max(abs(W),na.rm=TRUE)
if (is.null(breaks)) breaks <- round(seq(0,Wx,length=31),2)
if (fast)
map(W,xlim=range(lon(W)),ylim=range(lat(W)),
colbar=colbar,new=FALSE)
else {
attr(W,'variable') <- NULL
attr(W,'unit') <- NULL
map(W,xlim=range(lon(W)),ylim=range(lat(W)),projection='sphere',
colbar=colbar,verbose=verbose)
}
figlab(paste(it,collapse='-'),ypos=0.99)
}Function that applies tests to theresults through applying PCA
pcaasses <- function(Y) {
## Cut away blocks with missing data:
nv <- apply(coredata(Y),2,nv)
Y <- subset(Y,is=nv > 12000)
Y4s <- as.4seasons(Y,nmin=60)
y1 <- subset(Y4s,it='Jan')
y2 <- subset(Y4s,it='Apr')
y3 <- subset(Y4s,it='Jul')
y4 <- subset(Y4s,it='Oct')
## Fill in missing
y1 <- pcafill(y1)
y2 <- pcafill(y2)
y3 <- pcafill(y3)
y4 <- pcafill(y4)
## Apply PCA:
pca.djf <- PCA(y1)
pca.mam <- PCA(y2)
pca.jja <- PCA(y3)
pca.son <- PCA(y4)
## Examine the results:
print('DJF - PCA')
plot(pca.djf,new=FALSE)
print('MAM - PCA')
plot(pca.mam,new=FALSE)
print('JJA - PCA')
plot(pca.jja,new=FALSE)
print('SON - PCA')
plot(pca.son,new=FALSE)
## compare the station data with reanalysis - based on CCA
## Read reanalysis:
X <- retrieve('air.mon.mean.nc',lon=range(lon(Y))+c(-3,3),lat=range(lat(Y))+c(-3,3))
X4s <- as.4seasons(X)
eof.djf <- EOF(subset(X4s,it='Jan'))
eof.mam <- EOF(subset(X4s,it='Apr'))
eof.jja <- EOF(subset(X4s,it='Jul'))
eof.son <- EOF(subset(X4s,it='Oct'))
cca.djf <- CCA(pca.djf,eof.djf)
cca.mam <- CCA(pca.mam,eof.mam)
cca.jja <- CCA(pca.jja,eof.jja)
cca.son <- CCA(pca.son,eof.son)
print('DJF - CCA')
plot(cca.djf,new=FALSE)
print('MAM - CCA')
plot(cca.mam,new=FALSE)
print('JJA - CCA')
plot(cca.jja,new=FALSE)
print('SON - CCA')
plot(cca.son,new=FALSE)
## The quick and effiecint way to downscale all the stations:
ds.djf <- DS(subset(pca.djf,pattern=1:4),eof.djf)
ds.mam <- DS(subset(pca.mam,pattern=1:4),eof.mam)
ds.jja <- DS(subset(pca.jja,pattern=1:4),eof.jja)
ds.son <- DS(subset(pca.son,pattern=1:4),eof.son)
print('DJF - DS')
plot(ds.djf,new=FALSE)
#dev.new()
print('MAM - DS')
plot(ds.mam,new=FALSE)
#dev.new()
print('JJA - DS')
plot(ds.jja,new=FALSE)
#dev.new()
print('SON - DS')
plot(ds.son,new=FALSE)
return(list(pca.djf=pca.djf,pca.mam=pca.mam,pca.jja=pca.jja,
pca.son=pca.son,y.djf=y1,y.mam=y2,y.jja=y3,y.son=y4))
}Maximum temperature
The code below shows how the functions are used and applied to the data for the daily maximum \(T_x\). The first chunck of code carries out some preliminary analysis, the second shows how the results are derived, and the third show how the plots are made.
#-----------------------------------------------------------------------
# Define season and parameter
param <- 'tmax'
FUN <- 'mean'
reanalysis <- 'air.mon.mean.nc'
#reanalysis <- 'ERAINT_t2m_mon.nc'
FUNX <- 'mean'
data.path <- 'dse.Tz.NCEP/'
predictands(param) -> Y
Y <- subset(Y,it=c(1961,2011))
lon=range(lon(Y))+c(-10,10); lat=range(lat(Y))+c(-15,15)
## Check the predictands:
X <- pcaasses(Y)## [1] "DJF - PCA"
## [1] "MAM - PCA"
## [1] "JJA - PCA"
## [1] "SON - PCA"
## [1] "Warning : Calendar attribute has not been found in the meta data and will be set automatically."
## [1] "DJF - CCA"
## [1] "MAM - CCA"
## [1] "JJA - CCA"
## [1] "SON - CCA"
##
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%[1] "DJF - DS"
## [1] "MAM - DS"
## [1] "JJA - DS"
## [1] "SON - DS"
The PCA suggest that most of the variance in the station data for daily maximum temperature \(T_x\) can be represented by a small number of leading modes. For the SON season, however, the second mode is more prominent, suggesting more complicated structure where spatial variability has less coherency between the stations.
The CCA suggests dissimilar spatial temperature pattern with the higest temporal correlation in the predictand and predictor. The closest spatial similarity is seen in SON.
The DS analysis applied to the reanalysis suggests strongest match in MAM and weakest in DJF. The spatial patterns differ to some extent and there is not a great match between the predictor and predictand.
## Get the predictand -> Y
for (it in c('djf','mam','jja','son')) {
## Get the large-scale predictor:
if (!exists('predictor')) {
T2M <- retrieve(reanalysis,lon=lon,lat=lat)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
} else if (length(month(subset(predictor,it=it)))==0)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
#print(paste('Generating dse.',param,'.tz.rcp45.',it,'.rda',sep=''))
## Carry out the downscaling:
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''))) {
#dev.new()
dse.t2m.tz.rcp45 <- downscale(Y,predictor,it,param,FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''),dse.t2m.tz.rcp45)
}
#print(paste('Generating dse.',param,'.tz.rcp26.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp26 <- downscale(Y,predictor,it,param,rcp='rcp26',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''),dse.t2m.tz.rcp26)
}
#print(paste('Generating dse.',param,'.tz.rcp85.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp85 <- downscale(Y,predictor,it,param,rcp='rcp85',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''),dse.t2m.tz.rcp85)
}
}## [1] "Warning : Calendar attribute has not been found in the meta data and will be set automatically."#print('--- Completed downscaling ---')The March-May season.
## analysis:
load(paste(data.path,'dse.',param,'.tz.rcp85.mam.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
plot(subset(z,is='Arusha'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4g.',param,'.pdf'))
plot(subset(z,is='Musoma'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4h.',param,'.pdf'))
gridmap(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.pdf',sep=''))
validate(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.diagnose.pdf',sep=''))The plume plots look OK for Arusha and Musoma, and the observations fit reasonably well within the envelope of downscaled ensemble (RCP8.5). The gridded map suggests strongest warming in the north (~5C increase in 2099).
The rank test (Wilcox) indicates reasonable scores, although the rank distribution is not entirely flat. The diagnostic (target plot) for the diagnosis suggest that the trend in the ensemble is higher than in observation (horizontal axis) and that the inter-annual variance (vertical axis) is smaller in the ESD results.
The December-February season
load(paste(data.path,'dse.',param,'.tz.rcp85.djf.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.pdf',sep=''))
validate(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.validate.pdf',sep=''))
diagnose(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.diagnose.pdf',sep=''))A weaker warming (~2C) is projected in DJF, also with strongest response in the north. But the rank scores are weak.
The June-August season
load(paste(data.path,'dse.',param,'.tz.rcp85.jja.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.pdf',sep=''))
validate(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.validate.pdf',sep=''))
diagnose(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.diagnose.pdf',sep=''))The assessment suggest the JJA season is better downscaled, with reasonable rank scores. The projected warming is more uniform over Tanzania, albeit with some local variations.
The September-November season
load(paste(data.path,'dse.',param,'.tz.rcp85.son.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.pdf',sep=''))
validate(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.validate.pdf',sep=''))
diagnose(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.diagnose.pdf',sep=''))The SON season is associated with strongest warming in the north, and displays a clear north-south gradient. The warming is projected to be 0.7–1.3C by 2099, and the rank scores are moderately good - the distribution of the ranks had a U-shape rather than being flat.
Check a random station: is is well-represented by the PCA?
## Check the PCA-results:
z.djf <- as.station(X$pca.djf)
y.djf <- X$y.djf
par(new=FALSE)
plot.zoo(subset(z.djf,is=1),subset(y.djf,is=1),
main=paste('Test: Winter temperature at',loc(subset(y.djf,is=1))))
plot(subset(z.djf,is=1),lwd=3,col='grey70')
lines(subset(y.djf,is=1))The single station data could be well-recovered from the PCA product.
Implementation - Minimum temperature
Get and prepare the predictand data for the daily minimum temperature \(T_n\).
#-----------------------------------------------------------------------
# Define season and parameter
param <- 'tmin'
predictands(param) -> Y
Y <- subset(Y,it=c(1961,2011))
lon=range(lon(Y))+c(-10,10); lat=range(lat(Y))+c(-15,15)
## Check the predictands:
X <- pcaasses(Y)## [1] "DJF - PCA"
## [1] "MAM - PCA"
## [1] "JJA - PCA"
## [1] "SON - PCA"
## [1] "Warning : Calendar attribute has not been found in the meta data and will be set automatically."
## [1] "DJF - CCA"
## [1] "MAM - CCA"
## [1] "JJA - CCA"
## [1] "SON - CCA"
##
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%
|
| | 0%
|
|================ | 25%
|
|================================ | 50%
|
|================================================= | 75%
|
|=================================================================| 100%[1] "DJF - DS"
## [1] "MAM - DS"
## [1] "JJA - DS"
## [1] "SON - DS"
Apply the ESD models to implement thedownscaling
## Get the predictand -> Y
for (it in c('djf','mam','jja','son')) {
## Get the large-scale predictor:
if (!exists('predictor')) {
T2M <- retrieve(reanalysis,lon=lon,lat=lat)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
} else if (length(month(subset(predictor,it=it)))==0)
predictor <- subset(as.4seasons(T2M,FUNX=FUNX),it=it)
print(paste('Generating dse.',param,'.tz.rcp45.',it,'.rda',sep=''))
## Carry out the downscaling:
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''))) {
#dev.new()
dse.t2m.tz.rcp45 <- downscale(Y,predictor,it,param,FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp45.',it,'.rda',sep=''),dse.t2m.tz.rcp45)
}
print(paste('Generating dse.',param,'.tz.rcp26.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp26 <- downscale(Y,predictor,it,param,rcp='rcp26',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp26.',it,'.rda',sep=''),dse.t2m.tz.rcp26)
}
print(paste('Generating dse.',param,'.tz.rcp85.',it,'.rda',sep=''))
if (!file.exists(paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''))) {
dse.t2m.tz.rcp85 <- downscale(Y,predictor,it,param,rcp='rcp85',
FUN=FUN,FUNX=FUNX,lon=lon,lat=lat)
save(file=paste(data.path,'dse.',param,'.tz.rcp85.',it,'.rda',sep=''),dse.t2m.tz.rcp85)
}
}## [1] "Generating dse.tmin.tz.rcp45.djf.rda"
## [1] "Generating dse.tmin.tz.rcp26.djf.rda"
## [1] "Generating dse.tmin.tz.rcp85.djf.rda"
## [1] "Generating dse.tmin.tz.rcp45.mam.rda"
## [1] "Generating dse.tmin.tz.rcp26.mam.rda"
## [1] "Generating dse.tmin.tz.rcp85.mam.rda"
## [1] "Generating dse.tmin.tz.rcp45.jja.rda"
## [1] "Generating dse.tmin.tz.rcp26.jja.rda"
## [1] "Generating dse.tmin.tz.rcp85.jja.rda"
## [1] "Generating dse.tmin.tz.rcp45.son.rda"
## [1] "Generating dse.tmin.tz.rcp26.son.rda"
## [1] "Generating dse.tmin.tz.rcp85.son.rda"print('--- Completed downscaling ---')## [1] "--- Completed downscaling ---"The PCA suggest that most of the variance in the station data for daily maximum temperature \(T_n\) can be represented by a small number of leading modes. For the SON season, however, the second mode is more prominent, suggesting more complicated structure where spatial variability has less coherency between the stations.
The CCA suggests dissimilar spatial temperature pattern with the higest temporal correlation in the predictand and predictor.
The DS analysis applied to the reanalysis suggests strongest match in MAM and weakest in DJF. The spatial patterns differ to some extent and there is not a great match between the predictor and predictand.
The March-May season
## analysis:
load(paste(data.path,'dse.',param,'.tz.rcp85.mam.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
plot(subset(z,is='Arusha'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4g.',param,'.pdf'))
plot(subset(z,is='Musoma'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('figure4h.',param,'.pdf'))
gridmap(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.pdf',sep=''))
validate(z)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('MAM')
#dev.copy2pdf(file=paste('map.',param,'.mam.diagnose.pdf',sep=''))The plume plots look good for Arusha and Musoma, and strongest warming (3.5-4C) is projected for the eastern part of Tanzania. The evaluation of the quality gave moderate results - some discrepancies in terms of the rank scores.
The December-February season
load(paste(data.path,'dse.',param,'.tz.rcp85.djf.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.pdf',sep=''))
validate(z)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('DJF')
#dev.copy2pdf(file=paste('map.',param,'.djf.diagnose.pdf',sep=''))Strongest warming projected for the southeast part of the country (~3C), and the evaluation suggest moderate scores. Again, still some discrepancies.
The June-August season
load(paste(data.path,'dse.',param,'.tz.rcp85.jja.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.pdf',sep=''))
validate(z)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('JJA')
#dev.copy2pdf(file=paste('map.',param,'.jja.diagnose.pdf',sep=''))Most pronounced warming is projected for the souther part (~3C), and the assessment indicates moderate quality with some discrepancies. The ensemble results tend to account for lower inter-annual variability than in the observations and weaker trend (except for one station with higher trend).
The September-November season
#load(paste(data.path,'dse.',param,'.tz.rcp85.son.rda',sep=''))
load(paste(data.path,'dse.',param,'.tz.rcp45.son.rda',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
gridmap(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.pdf',sep=''))
validate(z)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.validate.pdf',sep=''))
diagnose(z,new=FALSE)
figlab('SON')
#dev.copy2pdf(file=paste('map.',param,'.son.diagnose.pdf',sep=''))
z <- as.station(dse.t2m.tz.rcp85)
plot(subset(z,is='Arusha'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('SON')
plot(subset(z,is='Musoma'),target.show=FALSE,map.show=FALSE,legend=FALSE,xlim=c(1900,2100),new=FALSE)
figlab('SON')
Strongest projected warming in the south (~3C) and moderate quality.
Check a random station: is is well-represented by the PCA?
## Check the PCA-results:
z.djf <- as.station(X$pca.djf)
y.djf <- X$y.djf
plot.zoo(subset(z.djf,is=1),subset(y.djf,is=1),
main=paste('Test: Winter temperature at',loc(subset(y.djf,is=1))))
plot(subset(z.djf,is=1),lwd=3,col='grey70')
lines(subset(y.djf,is=1))Single station records can be recovered from the PCA.
Summary
There are problems with some stations, which is visible in the PCA, CCA and DS carried out with just reanalysis. The spatial patterns indicate non-smooth and incoherent structures, suggesting problems with some of the stations and a mismatch with the reanalysis. Furthermore, the cross-validation reveal moderate scores (r typically less than 0.7), and the evaluations based on trends, range of inter-annual variability and rank score reveal some discrepancies.
The NCEP/NCAR reanalysis is probably better than ERAINT for ESD model calibration.